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Abstract 

An analytical solution to the non-equilibrium Marshak diffusion prob- 
lem in a planar slab of finite thickness is presented. Analytic expressions 
for the radiation and material energy densities as a function of space and 
time are derived using the Laplace transform method by summing over 
the first few residues at the poles of the transcendental equation. Inte- 
grated energy densities and leakage currents are also obtained in analytical 
form. Results for a planar slab of any finite thickness can be generated us- 
ing the analytic expressions of these quantities unlike the previous works 
wherein numerical results were generated to a specified degree of accuracy 
for a semi-infinite medium with semi analytical solutions. The benchmark 
results obtained in this work can be used to validate and verify non equi- 
librium radiation diffusion computer codes. 
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of finite thickness; Laplace transform method 
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1 Introduction 



The time dependent non-equilibrium radiation transport equation is non lin- 
early coupled to the material energy equation. Also the material properties 
have complex dependence on the independent variables. As a result, the time 
dependent thermal radiation transport problems are commonly solved numeri- 
cally. Benchmark results for test problems are necessary to validate and verify 
the numerical codes. Analytical solutions producing explicit expressions for the 
radiation and material energy density, integrated densities, leakage fluxes, etc. 
are the most desirable. 
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In the past, considerable amount of efforts have been applied for solving 
the Radiation Transport problem analytically [I]. All the available results are 
semi-analytical and the desired quantities are tabulated to some fixed accuracy 
level for a particular value of material parameters. The first semi-analytical 
solution was obtained by Marshak in which he considered a semi infinite planar 
slab with radiation incident upon the surface Assuming that the radiation 
and material fields are in equilibrium, similarity solutions were applied and 
the problem reduced to a second order ordinary differential equation which 
was solved numerically [3j. Later, Pomraning considered the non equilibrium 
radiation diffusion problem with the assumption that the material specific heat 
is proportional to the cube of the temperature so that the equations become 
linear in the radiative intensity and the fourth power of the material temperature 
[4]. Using Laplace transform method, Pomraning derived the semi analytical 
solutions for the surface quantities, integral quantities, and the distributions 
of radiative energy and material temperature as functions of space and time 
for the case when the speed of light is treated as infinity. Subsequently, the 
semi analytical results for finite speed of light was obtained by Su and Olson 
[5] . Semi analytical results were also obtained for the transport problem with a 
point source applied for a finite time [B] , [7] . The non-grey benchmark results for 
the two temperature non equilibrium radiative transfer were also generated [5] . 
The normal-mode expansion techniquefS] and spherical harmonics method[TU] 
have been applied for solving the steady state radiation transport equation in 
a finite slab. Numerical results for reflectivity and transmittivity in a finite 
planar slab having specularly reflecting boundaries has been obtained using the 
travelling wave transformation [TT] . 

In this paper, we solve the non equilibrium Marshak diffusion problem for 
a plane slab of finite thickness with radiation incident upon one of its surfaces. 
The difference from the semi infinite slab is that radiation will leak from the 
other surface due to its finite thickness and vacuum boundary condition. Non 
equilibrium diffusion codes can be more easily validated and verified against 
these benchmark results as there is no need to take a slab of very large size for 
avoiding boundary effects. As all real life systems are finite in size, these re- 
sults are also more practical and experimental results are more easily modelled. 
For the semi infinite slab, because of the multiple valuedness of the functions 
obtained by Laplace transform, inverting them using the inverse Laplace trans- 
form required evaluation of contributions from all the branch cuts. This resulted 
in integrals which had to be computed numerically [5 . The oscillations in the 
integrand resulted in difficulty in their convergence. The advantage of solving 
the finite problem is that because of the single valuedness of the Laplace trans- 
formed functions, the inversion using inverse Laplace transform is very simple. 
The sum of the residues at the singularities (poles) give the required solution. 
Thus analytical expressions for all the quantities of interest can be obtained 
for any slab width and parameter value. These results will be very useful for 
benchmarking and checking the sensitivity of non equilibrium radiation diffusion 
codes. 

The remainder of the paper is organised as follows. In Section [2j the ana- 
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Figure 1 : Flux incident on the left surface of a slab of thickness z = I 

lytical solution for the finite planar slab is derived. In Section [3j the results for 
the no retardation case and for finite speed of light are presented and explained. 
Finally, conclusions are given in Section 01 



2 Analytical solution 

We consider a planar slab of finite thickness which is purely absorbing and ho- 
mogeneous occupying < z < I. The medium is at zero temperature initially. 
At time t=0, a time independent radiative flux (Fi nc ) is incident on the sur- 
face at z=0 as shown in Fig. [T] Neglecting hydrodynamic motion and heat 
conduction, the one group radiative transfer equation (RTE) in the diffusion 
approximation and the material energy balance equation (ME) arc 1 

dE(z,t) d c dE{z,t) 



dt 0z 1 3k(T) dz 



= c K (T)[aT 4 (z,t) - E(z,t)] (1) 



Cv{T) ^W^ = - aT 4 (z,t)] (2) 

where E(z,t) is the radiation energy density, T(z,t) is the material temperature, 
k(T) is the opacity (absorption cross section), c is the speed of light, a is the 
radiation constant, and C V (T) is the specific heat of the material. 

The Marshak boundary condition on the surface at x=0 is given by 

mt) ' ~ c Fmc (3) 

where Fi nc is the flux incident upon the surface z=0. 
And that at z = I is 

The initial conditions on these two equations are 

E(z,0) = T(z,0) = (5) 
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To remove the nonlincarity in the RTE (Eq.flT])) and ME (Eq.©), opacity K is 
assumed to be independent of temperature and specific heat C v is assumed to 
be proportional to the cube of the temperature, i.e., C v = aT 3 . The RTE and 
the ME are recast into the dimensionless form by introducing the dimensionless 
independent variables given by 

]— Aden 

x=v3kz,t=( )t (6) 

a 

and new dependent variables given by 

_ c E(z,t) _ c aT 4 (z,t) 
u(x,t) = (-)[— },v(x,t) = {-)[— ] (7) 

^ * inc " inc 

With these new variables, the RTE and ME take the dimensionless form 
du(x, t) d 2 u(x,r) 



dr dx 2 

dv(x, t) 



+ v(x, t) — u(x, t) (8) 
= u(x, t) — v(x, t) (9) 



with the initial conditions 



u{x,0) = (10) 
v(x,0)=0 (11) 

And the boundary conditions on the surfaces are 

■ft„ + -^-0 (1 3, 

where b = \f^nl and the parameter e is defined as 

_ 16cr _ 4a 
ca a 

To solve Eqs. ([SJ - (IT5|) . we introduce the Laplace transform according to 



(14) 



f(s) = / dre-"- f(r) (15) 
Jo 

to obtain 

d 2 u(x : s) 

esu(x, s) — -7T — = v(x, s) — u(x, s) (16) 

ox z 

sv(x, s) = u(x, s) — v(x, s) (17) 
B(M) + ^.„ (19) 
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The solutions of Eqs. (fTHl) - p^|) in s space are obtained as 



3sin\B(s)(b - x)] + 2y/3/3(s)cosW(s)(b - x)\ 

u = = 

s[3sin(p(s)b) + 4V3/3(s)cos(P(s)b) - AP{s) 2 sin{P{s)b)} 

3sin[P(s)(b - ap] + 2yflf3{s)cos[f3{s){b - x)] 

V ~ s(s + l)[3sin(P(s)b) + iV3/3(s)cos(P(s)b) - 4P(s) 2 sin((3(s)b)} 

where /3(s) is given by 



(20) 
(21) 



^(a) = -_[l+ e ( a + l)] (22) 

Before solving for the radiation and material energy densities by inverting u 
and v, we first obtain the small and large r limits of u(x, r) and v(x, r) from the 
large and small s limits of Eqs. (l20l) and (l2~Tj) respectively. Using the theorems 



we have 



Um s -, 0o [sf(s)] = lim T ^o[f(r)] (23) 
iim s _>o[s/(s)] = limr^oolfir)] (24) 



u(x, 0) = v(x, 0) = (25) 

. . , . 36 + 2\/3 — 3a; . . 

u(x,t — > oo) — > vix.T — > oo — > = — (26) 

36 + 4^3 

Thus according to Eq. (|25p . at the initial instant, both the material and radi- 
ation energy densities are zero inside the slab. Eq. (|2^|) asserts that at infinite 
time the radiation and material energy density equilibrate among themselves. 
However, because of the finite thickness of the slab, flux leaks out of the right 
edge so that the energy densities vary linearly along the length of the slab. 

The solutions for u(x,t) and v(x,t) follow from u(x,s) and v(x,s) by in- 
verting them using the Laplace inversion theorem 



/(t) = ^ / dse^f( S ) (27) 
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where the integration contour is a line parallel to the imaginary s axis to the 
right of all the singularities of f(s). Both u and v are single valued functions 
and hence there are no branch points. However, there are an infinite number of 
poles obtained from the roots of the transcendental equation 

3sin(p(s)b) + 4V3f3(s)cos(j3(s)b) - A0 1 {s)sin{p{s)b) = (28) 

or,tan(P(s)b) = (29) 

The roots of the transcendental equation can by obtained by graphical plotting 
(as shown in Fig. [5] for b=l). 
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Figure 2: Finding the roots of the transcendental equation tan((3(s)) = f(j3) = 

4/3 2 (s)-3 



The contour is closed in the left half plane so that the large semi circle gives 
a zero contribution. 

For ease of generating results, in this work, MATHEMATICA has been 
used to obtain the first 30 roots namely, 0, 1.22826, 3.61221, 6.54624, 9.60463, 
12.7025, 15.8174, 18.9409, 22.0696, 25.2014, 28.3354, 31.4709, 34.6076, 37.745, 
40.8831, 44.0216, 47.1606, 50.2999, 53.4395, etc [H]. 

Corresponding to each root of f3(s), there exists two values of s, i.e., two 
simple poles. The poles are obtained from solution of Eq. (|2"2"1) . The root /3(s) = 
gives poles at s=0 and s=-ll. The residue at s=-ll is zero and that at s=0 is 
redundant as we already have a simple pole at s=0 as seen in the denominator 
of Eq. (I2U1) . According to the residue theorem, J c dse ST f(s) = 2nix (sum of the 
residues at the singularities). Hence the residue at s=0 is 

3sin[l3( S )(b - x)] + 2V3f3( S )cos{(3( S )(b - x)j 

Hm s ^n e [s — 0) = 

s[3sin(J3(s)b) + 4 % /3^(s)cos(/3(s)6) - 4:f3 2 (s)sin(f3(s)b)] 

_ 36 + 2^3 - 3z 
3b + 4V3 

which has a dependence only on position x and no dependence on time r. Thus 
the asymptotic (steady state) solution for the radiation energy density is 

36 + 2V3-3x 

36 + 4V3 (30) 

du(x,r) , Ov{x,t) 



which is also obtained by equating q^ t ' and V ^ T > i n Eqs. © and © to 

2 u(x 



u( x t) 

zero, solving — ^%' T = and obtaining the values of the constants from the 
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BC given by Eqs. (IT2"j) and (fT3j) . The contribution to the time dependent part 
will come from the other poles. Adding all these contributions will give us the 
complete space and time dependence of the radiation energy density. Similarly, 
the residue at s=0 for the material energy density v is 

T . „. n , 3sin[P(s)(b - x)} + 2V3P(s)cos[p(b - x)} 

um s _)0 e (s — 0) = 

s(s + l)[3sin(j3(s)b) + Ay/3 p(s) cos (P(s)b) - 4f3 2 (s)sin((3(s)b)} 

= 35 +2V ~ 3 ~ 3X (31) 
3b + 4^/3 v ; 

and that at s=-l is 0. Hence, the asymptotic (steady state) solution for the 
material energy density is 

3b + 2V3-3x 
VM= 36 + 4V3 (32) 

which is the same as the radiation energy density. This shows that the radia- 
tion and material energy density attain equilibrium after sufficient time. The 
poles for the second root of the transcendental equation (3 = 1.22826 are s=- 
25.49448516 and -0.59174484. Now, the residue at a simple pole s = s n is 

Lim e sr {s _ s ) 3sin[p(s)(b - x)} + 2V3(3( s )cos[(3(s){b - x)] 
OT " S " 6 s[3sin(P(s)b) + 4^3/3 (s) cos (j3(s)&) - 4/3 2 (s)sm(/3(s)&)] 

As s — s n is a pole of w, both the numerator and the denominator tend to 
0. Applying L'Hospital's rule to both the numerator and the denominator, the 
residue at a pole s = s n is 

e s " T [3sm(/3(s„)(& - x)) + 2V3(3(s n )cos/3(s n )(b - x)] 
s n [(3b + 4V3 - 4p(s n ) 2 b)cos{/3(s n )b) - (4y/30(a n )b + 80(s n ))sin(0(s n )b)}^^ 

where 

dP(s n ) = e(s n + l) 2 + l m) 

ds 2( Sn + l)^[- s "^fa +1)) ] 1/2 

For each value of f3(s), there are two poles (i.e., s values), namely 

- — Y t • The solution is obtained by summing 

over the residues from the first few poles. 

, , 36 + 2^-32; 

U{X, Tj = 1= 

36 + 4^3 



E 



* T [3sin(p(s n )(b - x)) + 2v / 3/3(s„)cos(/3(s„)(6 - a?))] 



• n [(36 + 4V3- 4p(s n ) 2 b)cos(0{s n )b) - (4V3/3(s n )b + 8/3(s n ))sm(/3(s„)6)]^ 



ds 



(34) 
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Similarly, the solution for the material density is 

36 + 2^3 - 3x 



E 



v(x,t) - - 
36 + 4^ 

e s » T [3sm(^(s„)(6 - ar)) + 2V3^(s n )cos(3( Sn )(b - a:)] 



s n (s n + 1)[(36 + 4\/3 - 4/3(s ri ) 2 6)cos(/3( S „)6) - (4V3/3(s n )& + 8/9(« n ))«n(^(«„)6)]^=i 

(35) 



3 Results 

3.1 Solution for e=0 

We first consider the e=0 case which arises when the speed of light is taken to be 
infinite so that radiation is not retarded initially. At infinite time, the radiation 
and material energy densities assume the same spatial dependence as for £ ^ 
case. 

, . , . 36 + 2%/3-3x , . 

ulx.T — > oo — > vix.T — > oo) — > 1= — (3d 

V ' V ; 3b + 4^3 K ' 

However, for r = 0, as s — >■ oo for e = 0, we obtain f3 = i where i = \] — 1. Thus, 

3sinh(b — x) + 2\/3cosh(b — x) , . 

ufa;, 0) = = (37) 

7sinh(b) + 4:V3cosh(b) 

v(x,0)=0 (38) 

Thus the material energy density is zero at r = as predicted by the initial 
condition. However, because of the absence of retardation effects, the radiation 
energy density attains a finite value consistent with the incoming flux of radi- 
ation. This behaviour is in agreement with that obtained in the case of a semi 
infinite planar slab for the no retardation case. 

We now obtain the solution u(x,t) and v(x,t) by inverting the Eqs. (|20l) 
and (|21l) using inverse Laplace transform as in the general case with e = 0. The 
difference from the e ^ case is that only one pole is obtained corresponding 
to a value of beta i.e., s = — gr^^ ■ in all the results generated in this paper, 
sum of the residues at the first 30 roots have been considered. Also, the value 
of b is chosen to be 1. 

The material energy density increases with time from a zero initial value 
whereas the radiation energy density attains a finite value even at very early 
times due to the absence of retardation effects as seen from Figs. [3] and 2] Both 
the energy densities attain the steady state value at r = 10 and maintains that 
value at a later time r = 100. 
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Scaled distance (x) 



Figure 3: Radiation energy density vs position in the slab at different times for 
e = 




Scaled distance (x) 



Figure 4: Material energy density vs position in the slab at different times for 
e = 
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Figure 5: Radiation energy density vs position in the slab at different times for 
e = 0.1 

3.2 Solution for e ^0 

For a finite value of e, both the radiation and material energy density increase 
gradually from zero as shown in Figs. [5] and [6] (for e = 0.1) . 

The first derivatives w.r.t. position of the radiation and material energy 
density are obtained as 



du(x, t) 



-3 



dx 36 + 4V3 



e s " T l-3l3(s n )cos(P( Sn )(b - x)) + 2^ 2 { Sn )sin{p{s n ){b - x))] 
^ s„[(36 + 4V3 - 4/3 2 ( s „)&)cos(/3(s»)6) - (4^/?(s„)6 + 8/3(s n ))sm(/?(s„)&)]^^ 



(39) 



and 



9w(x, r) 



9x 36 + 4\/3 



£ 



e s " T [-3/3( Sw )cos(/3(s„)(6 - a;)) + 2^/3 2 (s ra )sm(/3(g ra )(6 - x))] 

•> n (s n + 1)[(36 + 4^ - 4/3 2 (s n )6)cos(/3(s„)6) - (4V3/3(s n )6 + 8p(s n ))sin(p(s n )b)} 



and plotted in Figs. [7] and HI The gradient of both radiation and material 
energy densities obtain a constant value of — = —0.30217 after infinite 

OJ 3+4V3 
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Figure 6: Material energy density vs position in the slab at different times for 
e = 0.1 



time showing that there is a constant leakage of flux from the right surface due 
to the finite thickness. This result is different from the semi-infinite slab result 
where at infinite time, the entire halfspace is at a constant temperature with a 
uniform radiation field and hence there is no gradient and no flux [4]. 

The current of radiation leaking out from the left and right surfaces of the 
slab are 

The leakage currents are plotted as a function of time in Fig. [S] At later times, 
the values are seen to converge to the asymptotic values of 0.30217 and 0.6978 
respectively from the left and right surfaces. 

The averaged or integrated radiation energy density is given by 

f\ f b 3b + 2V3_3x 
ip r [T) = I u(x,T)dx = / = — ax 



36 + 4^ 



i [(36 + 4V3-4/? 2 ( s „)&) COS (/3(s„)6) - (4V|9(«„)& + 8 J 8(a n ))«n(0(«„)6)]^ 



ds 



b 



x / [3sin(i3(s n )(b - x)) + 2V30(s n )cos(/3(s n )(b - x))]dx 



o 
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Figure 7: Space derivative of radiation energy density vs position in the slab at 
different times for e = 0.1 
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Figure 8: Space derivative of material energy density vs position in the slab at 
different times for e = 0.1 
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Scaled time (t) 



Figure 9: Leakage currents from the two surfaces of the slab 



36 + 2\/3 & 3 b 2 

36 + 4V3 36 + AV3 2 



4(36 + 4^3 - 4f3 2 {s n )b)cos(P(s n )b) - (4V3/3( S „)6 + 8/3(s n ))sin(p(s n )b)] ^-1 



x [^(1 - cos{p{s n )b)) + 2VSsin(0(s n )b)] 

P( s n) 

Similarly the averaged or integrated material energy density is 

rb 



^m{r) = I v(x,r)dx 
Jo 



(43) 



36 + 2V3 b 3 b 2 

36 + 4^3 36 + 4^/3 2 



E 



t 1 *n{sn + 1)[(36 + 4V3 - 4/3 2 ( s „)6)cos(/3(s„)6) - (4>/30(s»)& + 8/3(s„))sm(/3(s„)6)]^f 



x i^T^i 1 - cos(p(s n )b)) + 2v / 3sm(/3( S „)6)] 
For b=l, the steady state integrated value is 0.5 as seen from FigfTUl 



(44) 
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10 



Scaled time (t) 



Figure 10: Integrated radiation and material energy densities as a function of 
time 



To check the consistency of the final results, we add Eqs. 
integrate over x from to b, yielding 



and 



and 



b 8u(x,t) 



dv(x, t) 
dr 



)dx 



d 2 u(x 1 r) 
dx 2 



dx 



du(b,r) du(0,r) 



dx 



dx 



i.e., 



d^ r (r) 9-0 m (r) du(b,r) <9«(0,t) 



Ot 



Ot 



dx 



dx 



(45) 



Using the expressions for the energy densities, their first derivatives in space and 
the integrated quantities, we find that both the left and right hand sides reduce 



toE„ 



[(36+4 % /3-4/3 2 (s„)f))cos(/3(s„)&)-(4\/3^(s„)6+8/3(s„))sm(0(s„)6)]^|^ " L pC^> 

cos(/3(s n )b))+2y/3sin(/3(s n )b)](e+ s \ 1 ) proving the consistency of the obtained 
solutions. 

As there are infinite number of residues, the exact solution is obtained only 
on adding all of them. However, the contribution from the poles decrease very 
sharply. To study convergence, we plot percentage error as a function of number 
of roots of the transcendental equation considered. As seen from Fig. [TTJ 2.1 
% error in the value of u(0, 2.5) is observed on considering only the first two 
roots i.e., the steady state result and residue for the two non zero poles. The 
errors arising due to non inclusion of higher order terms is more initially as the 
higher order poles contribute only at very small times because of the exponential 
term. The error falls sharply to a negligible value (0.005%) on considering the 
contribution from the first 6 roots i.e., first 11 poles. More accurate results can 
be obtained by adding residues from higher order poles. 
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Figure 11: Percentage error in the radiation energy density as a function of 
number of roots considered (N) 

4 Conclusions 

In this paper we have derived completely analytical expressions for the space 
and time dependent radiation and material energy densities inside a planar 
slab of finite thickness under the diffusion approximation. All other quantities 
of importance like the integrated values and the leakage currents are obtained 
directly from these analytical expressions for any parameter value like slab thick- 
ness, specific heat, etc. The series solutions are found to converge quickly and 
depending on the required degree of accuracy, the number of poles to be consid- 
ered is decided. These results can serve as benchmark for time dependent non 
equilibrium one dimensional radiation diffusion codes. 
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